#!/bin/bash

## clean.fq.gz -> sorted.markdup.bam

#Shell options
working_dir=$1
sample_list=$2
n_jobs=$3

fq_to_bam(){
#Tools
trimmomatic=/picb/lilab/tools/Trimmomatic-0.39/trimmomatic-0.39.jar
adapters=/picb/lilab/tools/Trimmomatic-0.39/adapters/
bwa=/picb/lilab/tools/bwa-0.7.17/bwa
samtools=/picb/tools/samtools-1.9/
gatk=/picb/lilab/tools/gatk-4.2.6.1/gatk

#Reference
reference=/picb/lilab5/share_data/ref/hg38/
bundle=/home/lilabguest7/lilab5_cxf/ref/gatk_bundle/ensembl_hg38

#Directory
if [ ! -d $2/cleanfq/ ]
then mkdir -p $2/cleanfq
fi 

if [ ! -d $2/info/ ]
then mkdir -p $2/info
fi 

if [ ! -d $2/qc/fastq_qc ]
then mkdir -p $2/qc/fastq_qc
fi 

## Filter adapters Trimmomatic
time java -jar ${trimmomatic} PE \
    $2/rawfq/$1_1.fq.gz $2/rawfq/$1_2.fq.gz \
    $2/cleanfq/$1.paired.1.fq.gz $1.unpaired.1.fq.gz \
    $2/cleanfq/$1.paired.2.fq.gz $1.unpaired.2.fq.gz \
    ILLUMINACLIP:$adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 

#Mapping bwa-mem
time ${bwa} mem -t 6 -M -Y -R "@RG\tID:$1\tSM:$1\tLB:$1\tPL:Illumina" \
             $reference/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
             $2/cleanfq/$1.paired.1.fq.gz $2/cleanfq/$1.paired.2.fq.gz | samtools view -Sb - > $2/bam/$1.bam 

#Sort samtools 
time $samtools sort -@ 4 -m 4G -O bam -o $2/bam/sort/$1.sorted.bam $2/bam/$1.bam 

#Mark Duplicates GATK
time $gatk MarkDuplicates \
    -I $2/bam/sort/$1.sorted.bam \
    -M $2/bam/markdup/$1.markdup_metrics.txt \
    -O $2/bam/markdup/$1.sorted.markdup.bam 

#Index samtools
time $samtools index $2/bam/markdup/$1.sorted.markdup.bam 

#BQSR
time $gatk BaseRecalibrator \
    -R $reference/Homo_sapiens.GRCh38.dna.primary_assembly \
    -I $2/bam/markdup/$1.sorted.markdup.bam \
    --known-sites ${bundle}/Mills_and_1000G_gold_standard.indels.hg38.ensembl.vcf.gz \
    --known-sites ${bundle}/dbsnp_146.hg38.ensembl.vcf.gz 

time $gatk ApplyBQSR \
    --bqsr-recal-file $2/bam/bqsr/$1.sorted.markdup.recal_data.table \
    -R $reference/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -I $2/bam/markdup/$1.sorted.markdup.bam 

#Index samtools
time $samtools index $2/bam/$1.sorted.markdup.BQSR.bam 
}

export -f fq_to_bam

parallel --jobs ${n_jobs} fq_to_bam :::: ${sample_list} ::: ${working_dir} 

